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The results of numerical simulations aimed at assessing the efficacy of Generalized 
Predictive Control (GPC) for active gust load alleviation using trailing- and leading- 
edge control surfaces are presented. GPC is a linear, time-invariant, multi-input/multi- 
output predictive control method that uses an AutoRegressive with exogenous input 
(ARX) model to characterize the system and to design the controller. The equations 
underlying the method are presented and discussed, including system identification, 
calculation of control law matrices, and calculation of the commands applied to the 
control effectors. Both embedded and explicit feedforward paths for incorporation of 
disturbance effects are addressed. Representative results from two types of simulations 
are presented. The first used a three-degree-of-freedom mathematical model of a mass- 
spring-das hpot system subject to a variety of user-defined external disturbances. The 
second used open-loop data obtained during a wind-tunnel test in which a wing model 
was subjected to sinusoidal vertical gust excitation; closed-loop behavior was simulated 
in post-test calculations. Results obtained from these simulations have been decidedly 
positive. In particular, the closed-loop simulations for the wing model showed reductions 
in wing root moments by factors as high as 1000, depending on whether the excitation is 
from a constant- or swept-frequency gust and on the direction of the response. 


Nomenclature 

AutoRegressive 

AutoRegressive with exogenous input 

disturbance acting on plant 

known and unknown portions of disturbance 

Generalized Predictive Control 

prediction and control horizons 

integer denoting time index 

number of data points used for system identification 

number of outputs 

order of AR disturbance model 

Observer Markov Parameters 

order of ARX model used for system identification 
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total number of inputs (r = r c + rj) 
number of control inputs 
number of disturbance inputs 
System IDentification 

coefficient matrices in multi-step output prediction equation 
input 

control input 

control input at time step k 

vectors of inputs, outputs, and disturbances at k lb time step 

vectors of past inputs, outputs, and disturbances at (k-j) th time step; j < k 

vectors of future inputs, outputs, and disturbances at (k+j) th time step; j > 0 

output 

matrix of observer Markov parameters 
control law gain matrices 
observer Markov parameters 


Introduction 

The Aeroelasticity Branch and the Army Research Laboratory’s Vehicle Technology 
Directorate at NASA Langley Research Center, in collaboration with Bell Helicopter Textron 
Inc, have recently completed a series of experimental evaluations of a predictive control 
technique known as Generalized Predictive Control (GPC) (ref. 1) to assess its potential for 
actively controlling the swashplate of tiltrotor aircraft to enhance aeroelastic stability in both 
helicopter and airplane modes of flight. The studies employed a 1/5-scale, semispan aeroelastic 
model of the V-22 that had been refurbished to serve as a tiltrotor research testbed for use in the 
Langley Transonic Dynamics Tunnel (TDT). The results of the exploratory experimental 
investigations conducted with this model using a multi-input/multi-output GPC-based method for 
active control are reported in references 2-5. The GPC algorithm employed was highly effective 
in increasing the stability of the critical (least stable) mode for all of the configurations and 
conditions tested. GPC has apparently not yet been employed to evaluate its potential for active 
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gust load alleviation (GLA) of fixed-wing aircraft using wing control surfaces. To this end, the 
intention was to upgrade and extend the GPC-based software system that was developed for the 
tiltrotor model and to conduct an investigation of the GPC method for active gust load alleviation 
as an adjunct to a planned test in the TDT of a large semi-span wing model of a SensorCraft 
concept being studied by Northrop Grumman as part of a High Lift over Drag Active (HiLDA) 
Wing Program (ref. 6). Although it was not possible to make all the necessary changes to the 
software in time for the subject test, appropriate data were taken during the test with the 
objective of conducting a post-test assessment and evaluation of the method. 

The purpose of this report is to present results of numerical simulations obtained from 
recently-developed MATLAB (ref. 7) programs intended to serve as the basis for the planned 
extensions to the tiltrotor GPC software system. First, a qualitative discussion of the essential 
features of GPC is given. Then, the equations underlying the method are presented and 
discussed, including system identification, calculation of control law matrices, and calculation of 
the commands applied to the control effectors. Procedures for incorporation of both known and 
unknown disturbance effects are addressed. Some general but relevant comments on the 
treatment of external disturbances, signal processing, and the codes developed and used for the 
simulations are then presented. Representative results from two types of simulations of open- 
and closed-loop behavior are presented next. The first set of simulations was made using a 
three-degree-of-freedom mathematical model of a mass-spring-dashpot system subject to a 
variety of user-defined external disturbances. The second set of simulations was made using 
data obtained during the HiLDA wing test conducted in the TDT in September 2004 (ref. 6) in 
which the model was subjected to sinusoidal vertical gust excitation at constant and continuously 
varying frequencies. The particular results presented from each type of simulation were chosen 
to illustrate the potential gust-load-alleviation capability that can be obtained with a GPC-based 
active control system. 


Basic Features of GPC 

Predictive control refers to a strategy wherein the decision for the current control action is 
based on minimization of a quadratic objective function that involves a prediction of the system 
response at some number of time steps into the future. A variety of predictive controllers have 
been proposed (ref. 8). Among these, Generalized Predictive Control, which was introduced in 
1987 (ref. 9), has received notable attention by researchers. GPC is a linear, time-invariant, 
multi-input/multi-output predictive control method that uses an AutoRegressive with exogenous 
input (ARX) model (i.e., a linear difference equation) to describe the input-output relationship of 
the system. The coefficient matrices of the ARX equation are determined by system 
identification. Such a model gives the current outputs as a linear combination of past inputs and 
outputs. The input-output equation is used to form a multi-step output prediction equation over a 
finite prediction horizon. An expression for the control to be imposed at the next time step is 
determined by minimizing the deviation of the predicted controlled plant outputs from the 
desired (or target) outputs, subject to a penalty on control effort. A version of the GPC 
procedure was developed at NASA Langley Research Center in 1997 for efficient computation 
and unknown disturbance rejection by Dr. Jer-Nan Juang and his associates. Their work resulted 
in a suite of MATLAB m-files that were collected into a predictive control toolbox that can be 
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used by researchers for GPC studies. A summary of the theory underlying their development is 
found in references 1 and 10-20, among others. 

The essential features of a GPC-based control system are depicted in figure 1 . The system 
has r c control inputs u, m measured outputs y, and is subject to external disturbances d that may 
be regarded as having known (<4) and unknown (d u ) portions. Measurement noise (not shown) is 
also present. There are two fundamental steps involved in GPC implementation: (1) 
identification of the system, and (2) use of the identified model to design a controller. Any 
system identification technique that returns an ARX model for the system may be used. System 
identification is done on-line, but not necessarily at every time point, in the presence of the 
disturbances acting on the system. If all the external disturbances are unknown (or treated as 
unknown), the effects of the disturbance-induced output are automatically embedded in the 
identified model and there is no need for explicit inclusion of disturbance terms in the ARX 
model (ref. 1). This approach leads to a feedback controller with embedded (or implicit) 
feedforward. However, if any of the disturbances acting on the system are either known or 
measurable, that data can be used in a feedforward path simultaneously with the feedback data to 
enhance closed-loop performance, as suggested in figure 1 . In this situation, the identified ARX 
model includes the effects of the disturbance-induced output associated with the unknown 
disturbances while an AutoRegressive (AR) model identified on-line can be used to include the 
effects of the measurable disturbances. This approach leads to a feedback controller with both 
embedded and explicit feedforward. In either case, the parameters of the identified model are 
used to compute the predictive control law. A random excitation um is applied initially with u c 
equal to zero to identify the open-loop system. Dither is added to the closed-loop control input 
u c if it is necessary to re-identify the system while operating in the closed-loop mode. It should 
be pointed out that this method lends itself to implementation strategies that allow all 
computations to be done on-line and in real time, as was demonstrated in experimental 
investigations with the tiltrotor model. 

The GPC software system developed for the tiltrotor model identifies the plant and 
computes the control law gain matrices assuming that the aerodynamic forces (the disturbances) 
acting on the model are unknown. This is the case of feedback with embedded feedforward 
described above and represents the capability that was incorporated in the m-files of the 
aforementioned predictive control toolbox. A Langley team has recently identified, implemented 
(as MATLAB codes), and verified the changes needed to include a feedforward path in the GPC- 
based software system developed for the tiltrotor model. This extension was motivated by the 
fact that the data needed for a feedforward path would be available from the HiLDA wing test in 
the form of gust-induced flow-angle measurements made upstream of the model. It is not 
difficult to envision other active control situations in which one or more of the disturbances 
acting on a system are measurable (or known) and which could take advantage of feedforward 
control. 


Pertinent Equations Underlying the GPC Method 

The key ingredients of the GPC method are system identification and output prediction 
based on the identified model; the control law is then derived by minimizing the error between 
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the predicted output and the desired output. The section begins by showing the general form of 
the system and control law equations. 

Form of System and Control Law Equations 

The relationship between the input and output discrete time histories of a linear, time 
invariant, multi-input/multi-output (MIMO) system is described by an ARX model that has the 
form 


y(k) - a l y{k-\) + a 2 y(k-2) + ••• +a p y{k-p) + /3 0 u(k) +/3 l u(k- 1 ) + ••• + J3 p u(k-p) 

+ y 0 d ( k ) + Y\d{k- 1 ) + • ■ ■ + y p d ( k-p ) 

This equation states that the current output y(k) at time step k may be estimated by using p sets of 
the previous output, input, and disturbance measurements, y(k-l), y(k-p), u(k-l), ...,u(k-p), and 
d(k-l),...d(k-p); and the current input and disturbance measurements u(k) and d(k). The integer p 
is called the order of the ARX model. The coefficient matrices a„ /f and 7 , are referred to as 
observer Markov parameters (OMP) or ARX parameters and are the quantities to be determined 
by the identification algorithm. Closed-loop performance is enhanced by performing the system 
identification in the presence of the external disturbances acting on the system, thereby ensuring 
that disturbance information will be incorporated in the system model. The goal of system 
identification is to determine the OMP based on input and output data. The OMP may be 
determined by any identification technique that returns an ARX model of the system. 

The ARX model is used to design the controller that, in the case of a regulator problem, 
has the general form given by 


u c (k) =a c x y{k- 1) + a^y(k- 2 ) + ••• +a p y(k-p) +ffiu(k- 1) +p c 2 u(k-2 ) + ••• +J3 p u(k-p) 

+ S L X d(k-l) +S c 2 d(k- 2 ) + • ■ • + S p d(k-p) +y c l d(k) +y c 2 d(k + 1 ) + • • ■ + y c p d(k + p - 1 ) * ” * 

Equation 2 indicates that the current control input u c (k) may be computed using p sets of the 
previous input and output measurements, u(k-l),..., u(k-p), d(k-l),..., d(k-p), and y(k-I ), .... y(k- 
p); and the future disturbances d(k), d(k+l), d{k+p-l). The coefficient matrices a, [f, d c \ and 

y c are the control law gain matrices. If the p past disturbances d{k-T), d{k-T),..., d(k-p) are 
measurable (known) they can be inserted directly into equation 2. The future disturbances, d(k), 
d(k+l), d(k+p-I), are typically unknown. However, if there is correlation in the disturbance 

signal (such as in the sinusoidal gust excitation of the HiLDA wing model) an AutoRegressive 
(AR) model of the disturbance can be determined and used to predict the future disturbances 
needed in equation 2. 

System Identification 

System identification (SID) in the presence of the operational disturbances is the first of 
the two major computational steps of the GPC method. The number of control inputs is r c , the 
number of disturbance inputs is r ( j, the number of measured outputs is m, and r = r c + rj. To 
initiate the identification process, the system is excited with band-limited white noise. These 
independent random excitations are applied to the r c inputs simultaneously and the m responses 
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are measured. The resulting input and output time histories ( u and y) at l time points are then 
used to fonn the data matrices y and V in the equation 

y = YV (3) 


where 


y= 


v(0) v(l) v( 2) • • -y(p) ■■■y(C-\) 


mxl 


( 4 ) 


and 


d( 0) d( 1) d( 2) • 

• • d(p) • 

• • d(£-\) 

m(0) m(1) m(2) 

u(p) 

u(l- 1) 

v(0) v(l) • 

.. v(p-V) ■ 

•• v(£-2) 

v(0) • 

•• v(p- 2) • 

■■ v(f-3) 


v(0) • 

•• v(£-p-'. 


[r + (r + m)p] x t 


Equations 4 and 5 follow from writing the equations for a discrete -time state-space observer 
model of a linear time-invariant system at a sequence of time steps k = 0, 1, 2, (f-1) and 

grouping (assembling) them into a matrix form (ref 11). Equations 4 and 5 are valid under the 
assumption that the system is initially at rest. For nonzero unknown initial conditions, the first p 
columns of y and V must be deleted to eliminate the effect of the initial conditions (ref. 11). The 
vector v(k) appearing in the data matrix V is formed from the vectors d(k), u(k) and y(k) 
according to 


v(k)= 


d(k) 
u(k ) 

y(k) 


(, r+m ) xl 


( 6 ) 


The order of the ARX model, p, and the number of time steps, t , are specified by the user. 
Some guidelines for their selection are given later. The sizes of key vectors and matrices are 
noted here and throughout the report as appropriate to facilitate computational considerations 
related to computer implementation. 


In forming the matrices given in equations 4 and 5, it has been assumed that a deadbeat 
observer has been added to the system (ref. 18). It is through this assumption that the matrix V 
can be reduced to a size amenable for practical numerical computation of its pseudo-inverse. 
The system identification process yields observer Markov parameters rather than system Markov 
parameters because of the inclusion of the observer. A complete discussion of these aspects of 
the development may be found in reference 1 1 . 

In equation 3, Y is the matrix of observer Markov parameters that is to be identified and 
has the form 
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( 7 ) 


Y= 


a. 


Y 0 Pq Y 1 P\ a \ Y 2 @2 a 2 P 3 ^3 ^3 

mx[r + (r + m) p] 


Yn P\ 


p 


The solution for Y is obtained by solving equation 3 for Y according to 

Y= yV'=yV r \w r ]-' < 8) 

where f denotes the pseudo-inverse. If the product V V T is a well-conditioned matrix of 
reasonable size, the ordinary inverse can be taken as shown. Otherwise, a pseudo-inverse must 
be used. 


Multi-Step Output Prediction Equation 

The one-step-ahead output prediction equation given in equation 1 is the starting point for 
deriving the multi-step output prediction equation that is needed for designing a GPC controller. 
Using equation 1, the output at time step k+j may be written in the form 

y(k + j) = a[ J) y{k-\) + a^y{k-2) + ■■■ +a ( J ) y{k-p) 

+ J3 0 u(k + j) +/3j ) u(k + j - 1) + • • • + jd{j ] u{k) +/3\ J) u(k-X) +/3j ) n(k - 2) + • • • +f3 ( p J) u(k-p) (9) 

+ y 0 d(k + j ) +y i 0 1) d(k + j - 1) + • • • + y ( 0 j) d(k) +y\ j) d{k-X) +y^d{k - 2) + • • • +y { f j ) d{k-p) 

where the coefficient matrices are given by recursive expressions involving the matrices a„ /?„ 
and y appearing in the ARX equation (ref. 10). The system identification process described 
earlier determines the matrices a h f:>„ and y Equation 9 shows that the output y(k+j ) at time step 
k+j may be estimated by using p sets of the previous output and input measurements, y(k-\ 
y(k-p ) and u(k-\), ...,u(k-p), and the (unknown) current and future inputs u(k), u(k+\), ...,u(k+j) 
and disturbances d(k), d(k+ 1 ), ...,d(k+j). The GPC algorithm is based on system output 
predictions over a finite horizon h p known as the prediction horizon. To predict future plant 
outputs, some assumption needs to be made about future control inputs. In determining the 
future control inputs for GPC, it is assumed that control is applied over a finite horizon h c known 
as the control horizon that is equal to or less than the prediction horizon. Beyond the control 
horizon the control input is assumed to be zero. Letting j in equation 9 range over the set of 
values j = 1, 2, ..., h p - 1, the resulting equations can be assembled into a multi-step output 
prediction equation having the fonn 

Y h (k)= T a, (k) + B u p (k-p) + A y p (k-p) 

, p , h p mxh c r c , ^ h p mxpr c pr c x\ h p mxpm pmx\ 

n P mx fcC /I Q\ 

+ V p d p (k-p ) + dj(k) 

hpitiXpi-j pr^ x 1 hpinxhp h p x 1 

The coefficient matrices % A, B, V p and V/ are formed from combinations of the observer 

Markov parameters a h />,, and y The quantity Vh p (k) is the vector containing the future outputs, 
whereas Uhdk) is the vector containing the future control inputs yet to be determined. The 
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quantities u p (k-p), y p {k-p), and d p (k-p ) are vectors containing the previous p sets of control inputs, 
outputs, and disturbance inputs, respectively. The vector df(k) contains the future disturbance 
inputs, which are usually unknown. However, if the disturbance is a coherent signal, future 
values of the disturbance can be estimated using an AR model of the disturbance obtained in a 
separate identification computation using a portion of the measured disturbance time history. 

The expanded form of this multi-step output prediction equation for h c ^ h p is shown in equation 

11 . 
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The observer Markov parameters /?„ and y ; determined in system identification are used 
to form the first block row in each of the coefficient matrices % A. B, V p and Vf in equation 10. 

These terms are colored red in equation 1 1 . The terms in the remaining block rows are computed 
using the recursive relations indicated in the boxes above the coefficient matrices in equation 1 1 
(ref. 10). All tenns in equation 1 1 are known, except for the h c sets of future commands, the h p 
sets of predicted responses, and the h p sets of future disturbance inputs tl/ (k). The goal of the 
GPC control algorithm is to determine the set of future commands u(k), u(k+ 1 ), ..., u(k+h c - 1) 
that are required to achieve a desired predicted response y(k), y(k+ 1 ), ..., y(k+h p - 1 ). It should be 
noted that the system Markov parameters, which are commonly used as the basis for identifying 
discrete-time models for linear dynamical systems, form the first block column in the matrix 7) 

the remaining block columns are formed from subsets of these parameters. The system Markov 
parameters are the pulse response of a system and are unique for a given linear system. The 
discrete-time state-space matrices A, B , C, and D are embedded in the system Markov 
parameters. 


Derivation of Control Law 


The predictive control law is obtained by minimizing the deviation of the predicted 
controlled response (as computed from the multi-step output prediction equation) from a desired 
response over a prediction horizon h p . To this end, one first defines an error function that is the 
difference between the desired (target) response vr (k) and the predicted response yy, (k): 


£=y T (k)-y hp (k) 

= y T ( k ) - T \ ( k ) - Bu p (k-p)- A Vp (k-p)- Vdp ( k-p ) - V f d f (k) 


( 12 ) 


An objective function J quadratic in the error and the unknown future controls is then fonned: 

J = £ T Qe + ul Ru h (13) 

n c n c 

Two weighting matrices are included in the objective function: R (symmetric and positive 
definite) is used to weight the control effort and stabilize the closed-loop system; Q (symmetric 
and positive semi-definite) is used to weight the relative importance of the differences between 
the target and predicted responses. Typically, R and Q are assumed to be diagonal with R having 
the same value w c along its diagonal and Q having the same value w r along its diagonal. 
Minimizing J with respect to Uh c (k) and then solving for u/Ak) gives 


hyX 1 


as the control sequence to be applied to the system over the next h c time steps. However, only 
the first r c values (corresponding to the first future time step) 
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( 15 ) 


U c (k) = ~Z c y T (k) + P c u p (k - p) + a c y p (k - p) + S'dfk - p) + f d f (k) 

r c x 1 


are applied to the r c control inputs, the remainder is discarded, and a new control sequence is 
calculated at the next time step. Equation 15 is the desired control law equation. The target 
response vr (k) is zero for a regulator problem and non-zero for a tracking problem. The 
coefficient matrices yf a, ff, d‘\ and y c are the control law gain matrices and have sizes r c x mh p , 
r c x pm, r c x pr c , r c x pr d , and r c x r d h p , respectively. 

Selection of Control Model Parameters 

An expression for estimating the order of the ARX model that is to be used for SID is 
given by 


p > ceil I 


number of system states + number of disturbance states 


m 


(16) 


where ceil denotes rounding up the value of the quantity in parentheses to the next higher 
integer. The number of system states is typically chosen to be twice the number of significant 
structural modes; the number of disturbance states is set to twice the (estimated) number of 
frequencies in the disturbance; m is the number of output measurements. If measurement noise 
is of concern, the value of p so computed should be increased to allow for computational poles 
and zeros to improve system identification in the presence of noise. In practice, simply choosing 
p to be about five times the number of significant modes in the system is often adequate. The 
prediction and control horizons are set according to the relations 

h p >p h c <h p (17) 

By extending the prediction and control horizons to very large values, the GPC solution 
approaches that of the linear quadratic regulator (LQR). Thus, GPC approximates an optimal 
controller for large h p (ref. 9). 

Weighting matrices R and Q are used to weight the control effort and to weight the relative 
importance of the differences between the target and predicted responses, respectively. In 
practice, the relative values of the weighting matrices play the major role in the optimization 
process. Therefore, the weighting matrix Q is often set to an identity matrix. The control weight 
matrix R then becomes the only one needed to be tuned to produce an acceptable solution 
without going unstable. Reducing w c increases controller performance but may saturate the 
actuators or drive the control system unstable. 

Practical Issues 

Two issues need to be addressed before applying the GPC method to a real problem: the 
treatment of external disturbances and basic signal processing guidelines. 
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Treatment of External Disturbances 


If the external disturbances acting on a system are unknown, the best that can be done is to 
identify the system in the presence of these disturbances so as to include their effects in the 
identified OMP. The control law obtained using the OMP with the embedded disturbance 
information is then used to compute the control inputs. This is the procedure that was employed 
successfully in the GPC-based active controls testing of the tiltrotor model. 

If any of the external disturbances acting on the system are known or can be measured a 
feedforward loop can be employed in combination with the feedback loop. There are two 
feedforward terms (e.g., see equation 10): a tenn involving past values d p and a term involving 
future values df If only past values are available one simply discards the other term. If future 
values are also available from an AR prediction model then both terms are retained. If there is 
correlation in the disturbance signals and the character of those disturbances is not changing 
significantly with time an AR model can be determined using an initial portion of the disturbance 
data. This model can then be used to predict the future values needed in equation 15. If the 
disturbances are changing significantly with time (e.g., due to frequency or amplitude variations, 
changes in the pattern or type of disturbance, etc) but there is coherence in the signals, an AR 
model can still be determined. However, for the AR model to adapt to the changing conditions it 
would have to be updated (ideally at each time step) using some type of recursive least squares 
algorithm. If there is no coherence in the disturbance signals no predictions of future 
disturbances are possible and only the tenn involving past values of the disturbance can be 
retained. These aspects are addressed in Appendix A. 

Basic Signal Processing Guidelines 

Several signal-processing considerations must be taken into account when using either 
measured data or simulated data in GPC-based active controls work. The measured response 
time histories must be passed through a low-pass filter with a cut-off frequency f c chosen so that 
the maximum frequency of interest is about 75% of f c . Assuming that f c is equal to the Nyquist 
frequency fy, the required sampling frequency f$ should be at least twice this value of f N to 
prevent aliasing. However, if fs is made too large the low-frequency modes will be poorly 
identified due to a loss of frequency resolution. A sampling rate between 2 to 3 times f\> is 
generally sufficient. Once the sampling frequency has been selected, the minimum number of 
data points required for system identification is at least 5-10 cycles of the lowest frequency mode 
in the measured response time histories. Normalization of the input and output data that is used 
for identification on the maximum actual or expected values of the data is often helpful 
numerically. 

It should be noted that active controllers are usually designed using a reduced-order 
mathematical model that includes only the most important subset of modes of a system. The 
frequency band spanned by the retained system modes sets the effective bandwidth of the 
resulting reduced-order controller. Such controllers may destabilize higher-order modes of the 
system due to truncation and spill over effects. The modes having frequencies in the transition 
region between retained and truncated modes are particularly susceptible. While these aspects of 
controller design are dealt with only casually below (e.g., through the choice of p , sampling 
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frequency, decimation factor, frequency range of interest, filter cutoff frequency, and number of 
data points used for SID), it is recognized that such considerations often require substantial 
attention by controls engineers. 


Results and Discussion 

Results from two types of closed-loop simulations are presented. One uses a 
mathematical model of a three-degree-of-freedom mass-spring-dashpot system subject to user- 
specified disturbance force time histories. The other uses gust response data obtained from the 
wind-tunnel test of the HiLDA wing model. The intent of these simulations is to demonstrate the 
efficacy of a GPC-based active control system for gust load alleviation and the improvement in 
performance associated with including a feedforward path if any of the external disturbances are 
known, measurable, or predictable. The results presented below have been selected primarily 
with these objectives in mind. 

Comment on Codes Used for Simulations 

Two MATLAB programs were developed and used for the simulation studies that 
provided the numerical results for this report. One is an extension of the three-degree-of- 
freedom (3-DOF) code that was used extensively to develop and evaluate the computational 
procedures that underlie the GPC software system that was used for the active controls testing of 
the tiltrotor model. This program treats a simple system with known dynamic properties and 
user-defined external disturbances. It also lends itself to easy modification should the need arise 
to address either a what-if question or a computational issue. The other program is new and was 
written specifically to access and process the data sets that were obtained during the wind-tunnel 
test of the HiLDA wing model (ref. 6). For that reason, it is tailored to treat that system. 
Although the programs are considerably different if viewed from a coding perspective, they are 
based on the same theoretical and computational principles and are therefore equivalent. 

Numerical Simulations Using 3-DOF Math Model 

The 3-DOF mass-spring-dashpot system used in the 3-DOF simulations is depicted in 
figure 2 and the corresponding equations of motion are summarized in Appendix B. The model 
has three masses, three springs, and three dashpots. Disturbance forces udi, ud 2 , and ud 3 and 
control forces uci, UC 2 , and UC 3 act on masses mi, m 2 , and m 3 , respectively. The degrees of 
freedom xj, x?, and X 3 are the absolute displacements of the corresponding masses. For the 
simulations, the masses were set equal to 1.0 lb-sec /in, the spring rates to 1000 lb/in, and the 
damping coefficients to 0.01 lb/in/sec. This yields natural frequencies of 2.240, 6.276, and 9.069 
Hz and corresponding damping ratios of 7.0367e-005, 1.9716e-004, and 2.8491e-004 for the 
three natural modes of the system. Disturbances, control forces, and responses (displacements, 
velocities, or accelerations) can be specified and selected in any combination (see Appendix B). 
However, for the results to be shown here, controls were imposed at masses 1 and/or 3 (except 
for one case of three control inputs), disturbances were imposed at masses 1 and/or 2, and 
accelerations were measured at masses 1 and/or 2. Both steady-state (constant amplitude and 
frequency) and non-steady-state (varying amplitude and frequency) disturbances were imposed 
on the system. 
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The steady-state disturbances (figure 3) were chosen to be sinusoidal with a frequency of 
6.276 Hz, which coincides with the second natural frequency of the system. The expressions 
used to represent these disturbances in the simulations were udj = 2*cos(2*pi*f*t) and ltd? = 
sin(2*pi*f*t). The non-steady-state disturbances had variable and/or discontinuous amplitude 
and frequency time histories. The first of these was composed of four variable-amplitude, 
constant-frequency segments (fig. 4a); the other represented a disturbance with increasing 
frequency and decreasing amplitude (fig. 4b). 

The MATLAB function c2d was used to convert the continuous-time state-space models 
of the 3-DOF system (Appendix B) to discrete-time using At =0.05 as set below. The resulting 

discrete-time state-space models, along with their corresponding discretized disturbances, were 
then used in the MATLAB function dlsim to generate the open-loop time histories needed for 
SID. 


The maximum frequency of interest for the 3-DOF system is just over 9 Hz. Based on 
the signal processing guidelines mentioned earlier, the sampling frequency/, was set to 20 Hz 

(A/=0.05 sec) and the number of time steps t of input/output data was set to 300. The control 

horizon h c was set equal to the prediction horizon h p and h p was set equal to the order of the 
ARX model p. The value of p was set to 6, 8, or 10 depending on the control case being 
considered and the number of disturbances acting on the system. For the cases of feedback-only 
and feedback with feedforward p was set to 6. The use of feedback with embedded feedforward 
requires p to be increased by (at least) two for each of the dominant frequencies present in the 
external disturbances acting on the system. This expedient provides the additional degrees of 
freedom needed to ensure that the effects of the unknown periodic disturbances are embedded in 
the observer Markov parameters that are computed during system identification. Thus, p was set 
to 8 when using the single disturbances in figure 4 and to 10 when using the double disturbances 
in figure 3. Based on previous experience with the 3-DOF model, the response weights w r were 
set to 1.0 and the control weights w c were taken to be equal and set to 0.3. In practice, the 
response weights w r are usually set to unity and w c is varied to tune the controller to the nuances 
of the subject system. Tuning proceeds by initially setting the value of w c to a “conservative” 
value (e.g., 1, as was initially done for the 3-DOF system) and then tracking the value of the 
objective function J given in equation 13 as w c is incrementally reduced until either a minimum 
value of J is identified or the control system actuators become saturated. This process for 
selecting w c is analogous to the familiar procedure for selecting the controller gain in a single- 
input/single-output system. Previous experience with the 3-DOF model indicates that the 
selected values are all reasonable but not optimum. No attempt was made to “optimize” (i.e., 
tune) the parameters p, h c , h p , w r , and w c to the nuances of the system under the different 
disturbances to elicit the best (lowest) possible closed-loop responses. In the simulations, the 
loop was closed after 100 time steps (5 sec) when using steady-state disturbances and 20 time 
steps ( 1 sec) when using non-steady-state disturbances. 

Results are presented below for three types of control cases: (1) feedback only, in which 
no account is taken of any external disturbances acting on the system during SID and controller 
design; (2) feedback with embedded feedforward, in which all the disturbances acting on the 
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system are assumed to be unknown while perfonning the SID and controller design; and (3) 
feedback plus feedforward, in which all of the external disturbances acting on the system are 
assumed to be measurable and are used during the SID and design of the feedforward controller. 
A second-order AR disturbance model (rid = 2) (see Appendix A) was used for predicting the 
future disturbances needed by the feedforward controller. This portion of the feedforward 
controller was adaptive and used RLS with a forgetting factor (ref. 21) of 0.9. Results are 
presented for both steady-state (constant amplitude and frequency) and non-steady-state (varying 
amplitude and frequency) disturbances. For ease of discussion, the second control case will be 
referred to simply as ‘embedded feedforward’ and the third control case will be referred to as 
‘explicit feedforward’. 

While some of the results shown below indicate the effect on performance of increasing 
the number of control effectors, there is no discussion of trade-offs between the number of 
control inputs used and the costs associated with implementation or power requirements. These 
would be important considerations in any studies aimed at full-scale implementation but are 
beyond the scope of this report. 

Representative Results 

Steady-State Disturbances : Representative results obtained from simulations made using the 
steady-state disturbances shown in figure 3 acting on masses 1 and 2 ( udi and ltd? in fig. 2) are 
given in figures 5-9. Loop closure (active control on) was taken to occur after five seconds. 
Figures 5-7 show results for feedback only, embedded feedforward, and explicit feedforward, 
respectively, for the case in which control is applied at mass 3 (wc?) and feedback is the 
acceleration of mass 1 ( y ). In all three cases, the response was quickly reduced to very small 
values when the loop was closed and then maintained at those low values thereafter. The closed- 
loop response levels for embedded feedforward (fig. 6) are essentially the same as those for 
explicit feedforward (fig. 7), but both are smaller than those for the feedback-only case of figure 
5. However, the feedback-only case also produces a considerably reduced response level. This 
result indicates that ignoring the disturbances during SID may produce an acceptable ARX 
model for designing a GPC controller for a system. The control inputs for all three cases are 
nearly the same. 

The effect of applying control forces at two masses ( uci and ucj) rather than just at one 
mass (ucj) as in figure 7 is illustrated in figure 8 for explicit feedforward. Comparing these 
results to those of figure 7 indicates that the use of an additional control force reduced both the 
closed-loop response and the magnitudes of the required control forces below the values 
associated with a single control, as might be expected. 

Figure 9 shows results for explicit feedforward for the case of control inputs at masses 1 
and 3 (ucj and tied as in figure 8 but with acceleration outputs at both masses 1 and 2 ( y andx, ). 
Both responses are very quickly quashed after the loop is closed. There is a slight increase in the 
magnitudes of the required control forces because two responses are now being worked by the 
controller rather than just one as in figure 8. 
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Non-Steady-State Disturbances : Typical results obtained from simulations made using the 
piecewise varying frequency and amplitude disturbance shown in figure 4a acting on mass 2 
(ud 2 ) and taking the feedback to be the acceleration of mass 1 ( x t ) are given in figures 10-14 for 

cases involving one to three control inputs. Loop closure (active control on) was imposed after 
one second. All responses were reduced to varying extents depending on the number of control 
inputs. Figures 10 and 11 show results using embedded feedforward and explicit feedforward, 
respectively, for the case in which control is applied at mass 3 (UC3). The closed-loop response 
using explicit feedforward is slightly better than for embedded feedforward while the magnitudes 
of the control inputs are either slightly larger, smaller or about the same depending on the 
segment of the disturbance that is active. 

Figures 12 and 13 show results for embedded feedforward and explicit feedforward, 
respectively, for the case of control inputs at masses 1 and 3 (uci and UC3). The closed-loop 
response is only slightly better using explicit feedforward (fig. 13) than embedded feedforward 
(fig. 12) while the required control inputs are about the same. Increasing the number of control 
inputs from one to two produces lower responses (compare figure 12 to figure 10 and figure 13 
to figure 11). 

Figure 14 shows results using explicit feedforward for the case in which control inputs 
are at masses 1, 2 and 3 ( uci , UC2 and UC3). The closed-loop responses and controls in this case 
are improved relative to the two-control-input results shown in figure 13, as should be expected. 

Results from simulations made using the disturbance shown in figure 4b acting at mass 2 
(udi) and taking the feedback to be the acceleration of mass 1 ( y ) are given in figures 15-17. 

Loop closure was taken to occur after one second. Figures 15 and 16 show results for embedded 
feedforward and explicit feedforward, respectively, for the case of a control input at mass 3 (wcj). 
The response is reduced in both cases but is decidedly better using explicit feedforward with 
adaptive identification of the disturbance. Control inputs for the latter are slightly larger during 
the early (lower frequency) portion of the disturbance but lower during the latter (higher 
frequency) portion of the disturbance. Figure 17 shows results using explicit feedforward for the 
case of control inputs at masses 1 and 3 (uci and UC3). The closed-loop response in this case is 
lower than that for the case of one control in figure 16, and control inputs are slightly smaller 
during the early (lower frequency) portion of the disturbance but generally the same or larger 
during the latter (higher frequency) portion. 

Summary of Observations : Taken as a whole, the 3-DOF simulations indicate that GPC-based 
controllers are able to significantly reduce acceleration responses under a variety of imposed 
disturbance conditions and to do so very quickly after closing the loop. For the cases involving 
steady-state disturbances (figs. 5-9), the lightly-damped system was being excited at one of its 
natural frequencies and the open-loop responses are all diverging with time. Thus, the 
percentage reduction in response after closing the loop in these cases depends on when the loop 
is closed as the magnitude of the steady-state closed-loop response is essentially independent of 
loop closure time. For figures 5-8 (2 disturbances, 1 output, 1 control) the amplitude of the 
open-loop response at time of closure (5 sec) is about 70 in/sec . With respect to this value, the 
reduction in acceleration responses ranged from about 90% to essentially 100%. For the case of 
figure 9 (2 disturbances, 2 outputs, 2 controls) both x, and x 2 were effectively reduced to zero. 
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The results in figures 10-17 are for disturbances with varying frequency and amplitude and can’t 
be quantified with a single number because the percentage reduction in response varies with the 
frequency and amplitude of the disturbance. However, depending on the position in the time 
history, the reduction varies from about 50% to essentially 100%. Given that the maximum 
values of the disturbance forces vary from one to two pounds in magnitude the required steady- 
state control forces in figures 5-17, which range from about 0.2 lbs to 2 lbs, are considered 
reasonable. 

With respect to the effects on performance associated with varying the number of control 
inputs and response outputs, the 3-DOF simulation results indicate that increasing the number of 
controls r c for a fixed number of outputs in improves perfonnance (both responses and control 
inputs are lower). Conversely, increasing the number of outputs m for fixed r c increases the 
required control inputs slightly but still quells all the responses. If both m and r c are increased 
simultaneously, all the responses are reduced with only a slight increase in the required control 
inputs. As mentioned earlier, no attempt was made to tune any of the GPC parameters to the 
nuances of the system under the different disturbances to obtain the best (lowest) possible 
closed-loop responses. However, the results (not shown) of some limited parametric variations 
indicate that increasing the order p of the ARX model and/or reducing the magnitude of the 
control weight w c can improve the performance of the controller in all of the 3-DOF simulations 
presented above. 

Numerical Simulations Using Data from Wind-Tunnel Test of HiLDA Wing Model 

The HiLDA wing model, as mounted in the Langley Transonic Dynamics Tunnel (TDT), 
is shown in figure 18. The model has four trailing-edge control surfaces that span the lifting 
portion of the wing and one leading-edge control surface located near the wing tip. A simulated 
vertical gust field was applied to the test section flow in the fonn of a sinusoidal oscillation of 
the flow direction. This oscillating flow was generated by an airstream oscillator system 
consisting of a biplane arrangement of vanes located on the side walls of the tunnel about 70 feet 
upstream from the model (fig. 19). The gust-induced flow angle was measured by a vane-type 
flow direction transmitter that was mounted on the top of a seven-foot high post attached to the 
tunnel floor six feet ahead of the model (fig. 18). A schematic of the wing planform is depicted 
in figure 20, which shows the arrangement of the five control surfaces that were used as control 
effectors and the locations of the three strain gage sensors that provided the feedback response 
time histories for the simulations. The reader is referred to reference 6 for a description of the 
model. 

Test Conditions 

The data used in the simulations presented below were obtained during a run in the TDT 
heavy gas (R134a) test medium at a dynamic pressure of 51 lb/ft 2 . The corresponding Mach 
number and velocity were 0.25 and 138 ft/sec, respectively. Sinusoidal vertical gust excitations 
were induced in the test section by the vertical components of the tip vortices shed from the 
oscillating biplane vanes. The amplitude of oscillation of the airstream oscillator was set to ±6 
degrees. This amplitude was sufficient to produce gust angles as large as one degree at the 
location of the model at the lower airspeed/frequency combinations. The airstream oscillator 
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excitation signatures consisted of constant frequency dwells at the natural frequencies of the first 
and second vertical bending inodes (flv ~ 2 Hz, f2v ~ 8 Hz), and frequency sweeps from 0.5 Hz 
to 10 Hz in 30 seconds. For SID purposes random excitations of the five control surfaces were 
imposed simultaneously and consisted of uncorrelated, normally-distributed random numbers 
with mean zero and a variance chosen to limit the oscillation angles of the control surfaces to 
about ± 2 degrees. 

Test Procedures 

The SID data for the case in which the gust excitation frequency is held constant were 
obtained as follows: (1) Set tunnel dynamic pressure to the desired value; (2) Set the frequency 
of the airstream oscillator to its desired value; (3) Impose a random excitation on each of the five 
control surfaces; and (4) Record the resulting input and output time histories for 35 seconds. For 
the case in which the gust frequency is swept from a low value to a high value the procedure 
was: (1) Set tunnel dynamic pressure to the desired value; (2) Set the frequency of the airstream 
oscillator to its starting value of 0.5 Hz; (3) Impose a random excitation on each of the five 
control surfaces; and (4) Sweep the gust frequency continuously from 0.5 Hz to 10 Hz in 30 
seconds while recording the resulting input and output time histories. 

Simulation Conditions and Parameter Settings 

The original data were taken using a sampling rate of 500 Hz and recorded for 35 
seconds, yielding data records consisting of 17500 time samples. To exclude transients, the data 
outside time samples 1000 and 15000 were discarded. Since the wing modes to be controlled 
were less than 15 Hz, the truncated data were then decimated by a factor of 10 using the 
MATLAB command decimate, resulting in data representing a sampling rate of 50 Hz. Finally, 
the mean value of each measurement was subtracted from its time history. 

The conditioned input-output data were used to identify the plant and to develop the 
controller. SID was accomplished using the command arx_ud from the predictive control 
toolbox. The outputs from strain-gage sensors SBI2, SBY, and STI2 shown in figure 20 were 
selected as the quantities to be minimized collectively (m = 3). All five control inputs, CTL1, 
CTL2, CTL3, CTL4, and CTL5, were utilized simultaneously as control effectors (r c = 5). The 
feedforward signal was taken from the flow direction transmitter. 

System identification returned the ARX parameters of equation 1 . The control law gain 
matrices in equation 15 were then obtained by processing the ARX parameters according to the 
computational methods described earlier in this report. The ARX parameters were used for the 
time history simulation of both the open- and closed-loop plant, as well as for the design of both 
the embedded feedforward and explicit feedforward controllers. Controller design was 
accomplished using a modified version of the m-file msop from the predictive control toolbox. 
The same controller was used for embedded feedforward and explicit feedforward simulations, 
the difference being that in the former case the feedforward controller path was not utilized. For 
closed-loop simulations employing embedded feedforward, only the tenns involving the past 
control inputs u p and feedback outputs y p in equation 15 were used to compute the control inputs 
u c . For cases employing explicit feedforward, all four terms were used. In order to use the df 
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term of equation 15, some estimate must be made of future disturbance inputs given the history 
of past disturbance measurements. As discussed earlier, depending on the character of the 
disturbances acting on the system, either an adaptive or a non-adaptive AR model can be used to 
predict future disturbances. Here, an adaptive 6 th order AR model was chosen to predict the 
future disturbances. A recursive least squares procedure with a forgetting factor of 0.97 was 
used to update the AR model parameters every time step, thus allowing for enhanced tracking of 
the swept sine gust. 

A comment is in order here regarding some additional data processing that was necessary 
to allow a proper comparison of measured open-loop time histories with simulated closed-loop 
time histories. The desired open-loop time histories (quiescent controls during each of the three 
gust-excitation conditions) were not among the data sets collected during the test. The closest set 
of data was that taken for SID wherein the five control surfaces were randomly excited during 
gust excitation. Fortunately, it is possible to compute the desired open-loop responses using that 
data set. Essentially, one first performs an SID using as input the random excitations of the 
control surfaces and the gust vane measurements. This yields the OMP a h fl, and y ; of equation 
1 . One then performs a time -history simulation for y(k) using equation 1 but retaining only the 
a, and y, terms on the right-hand-side and using for d the same set of gust vane measurements 
used for the SID. The resulting time histories are the desired open-loop time histories. To 
maintain consistency, the same set of gust vane measurements is added to the plant response y at 
each time step during the closed-loop simulations. 

The system model identified as described above was used in both the open-loop and 
closed-loop simulations for the cases of both embedded feedforward and explicit feedforward 
controllers. The control inputs u c {k ) were set to zero for the open-loop simulation. For the 
closed-loop simulation the control signal computed using equation 15 was used. The control 
input for embedded feedforward used only the past outputs y p and past control inputs u p ; no use 
was made of the past and future feedforward terms (d p and dj). For the explicit feedforward 
controller the computed control signal u c (k ) was obtained using all four terms on the right-hand- 
side of equation 15: the past outputs, past control inputs, past measured feedforward values, and 
future predicted feedforward values. The same weights that were used for the 3-DOF 
simulations (w r = 1 and w c = 0.3) were tried for the HiLDA wing data and found to be acceptable. 
The value of p was estimated using equation 16 and set to 10. Its value was not varied with 
controller type as for the 3-DOF model. The controller was turned on five seconds into the 
simulations. 

Representative Results 

Figures 21 through 62 show simulated open- and closed-loop responses of wing root 
vertical bending (SBI2), fore-and-aft bending (SBY), and torsion (STI2) to sinusoidal vertical 
gust excitation held constant at the natural frequencies of the first and second vertical bending 
modes and while sweeping the frequency of the gust from 0.5 Hz to 10 Hz. Bending and torsion 
moment response time histories, power spectral densities of bending and torsion moment time 
histories, and control angle input time histories are shown for each of the three gust excitation 
conditions for two control cases: (1) Feedback with embedded feedforward; and (2) Feedback 
with explicit feedforward and using an adaptive AR model for predicting the future disturbances. 
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Feedback with Embedded Feedforward : Figures 21-27 indicate controller performance in 
reducing wing root moments due to gust excitation at the natural frequency of the first vertical 
bending mode (2 Hz). The SBI2 closed-loop time history (fig. 21) shows a dramatic reduction of 
that modal response almost immediately after the controller is turned on. The corresponding 
SBY and STI2 closed-loop time histories (figs. 22-23) are also significantly, and promptly, 
reduced. The PSD plots (figs. 24-26) indicate that the primary effect has been to reduce the 
response at the frequency of excitation (the first mode), as would be expected. However, the 
closed-loop spectra are not flat indicating that there is residual correlation in the sensor signals. 
It may not be possible to remove all correlation in regulated systems due to stability issues, the 
presence of unobservable and/or uncontrollable modes, or limits on actuator authority. Based on 
experience with GPC applied to other systems, it is known that tuning of the GPC parameters 
(horizons, control penalty, order) oftentimes leads to enhanced performance. However, no 
adjustments were made in these parameters to optimize performance in the results shown here or 
below. Figure 27 shows the commanded control deflection angles to be small once transients 
decay. In a physical control system, the overshoot occurring at loop closure would be 
suppressed by low-pass filtering the control signal and/or placing upper and lower bounds on the 
allowable control amplitudes in both the hardware and software. The roll-off occurring in the 
magnitudes of the PSD curves at frequencies above about 20 Hz is due to the filtering action of 
the function decimate that was used to decimate the original wind-tunnel data. 

Figures 28-34 summarize controller performance due to gust excitation of the wing at its 
second vertical bending natural frequency (8 Hz). The response time histories (figs. 28-30) show 
a dramatic and essentially immediate drop in amplitude once the control is applied. The PSD 
plots (figs. 31-33) show a relatively flat controlled spectrum. Again, the primary effect has been 
to reduce the response at the frequency of excitation (the second mode), as would be expected. 
In this case the response in that mode is virtually eliminated. Additionally, the SBY responses at 
5 Hz (the first fore-and-aft bending mode) and at 16 Hz are greatly reduced. Figure 34 shows the 
control effort to be small once transients decay. In practice, the overshoot occurring at 5 seconds 
would be suppressed by low pass filtering the control signal and/or placing upper and lower 
bounds on the allowable control amplitudes. 

Figures 35-41 show controller perfonnance for the case in which the gust frequency is 
swept from 0.5 Hz to 10 Hz. In comparison to the constant-frequency gust excitation cases of 
figures 21-34, controller performance is only fair. This is due to the fact that SID in the presence 
of a varying-frequency gust leads to an internal (i.e., embedded) disturbance model that is 
effectively an average model over the bandwidth of the gust. This leads to a controller that can 
not effectively reduce the responses at all the frequencies contained in the gust as in the constant- 
frequency gust cases. One possible solution to this problem is to use a GPC controller in which 
all the coefficient matrices x, a c , ff, S c , and y in equation 15 periodically adapt (see, for 
example, ref. 15). Figure 41 shows that the control deflection angles in this case are large 
compared to those required for a constant-frequency gust excitation. It is thought that this is 
because the responses to the disturbance input (gust vane measurement added to y) were larger 
than the responses to the random disturbances input by the control surfaces during SID. There 
was thus a greater correlation in the gust vane disturbance input-to-output path than in the 
control input-to-output path. A more accurate mapping from control input u to system output y 


19 



would produce a better model of the plant. This could be achieved by either increasing the 
amplitude of the random noise applied to the control inputs for SID or reducing the level of the 
disturbance (gust vane measurement added to y) applied to the system during SID. 

Feedback with Explicit Feedforward : Figures 42-48 show the performance of the explicit 
feedforward controller for the case in which the gust frequency is set to that of the first vertical 
bending mode of the model. Recursive least squares is employed to update the coefficients in 
the AR model used to estimate the future values of the gust disturbance df given the past 
disturbance values d p , as described in Appendix A. In comparison to the embedded feedforward 
controller (figs. 21-27), the explicit feedforward controller exhibits better performance as can be 
seen by comparing time histories, PSDs, and required control angles. The steady-state control 
angles are all small after the transients have decayed. The improved performance is due to the 
incorporation of the measured past disturbance values d p into the controller as well as the ability 
of the future component of the feedforward path to tune itself to the changing conditions at each 
time step, thereby providing improved estimates of future disturbance values df. 

Figures 49-55 are for the case of explicit feedforward with a gust frequency set to that of 
the second vertical bending mode of the model. The adaptive part of the controller is the matrix 
Df p in equation A4 that is used to predict the future disturbance values d/(k) needed in the last 
tenn in equation 15. Enhanced performance over the embedded feedforward case is observed as 
it was for the first vertical bending mode. The control angles are again all small once the 
transients decay. It is of interest to note that the closed-loop spectra (figs. 53-55) are effectively 
flat, indicating that the controller is approximately optimal. 

Figures 56-62 show performance of the explicit feedforward controller for the case in 
which the gust frequency is swept from 0.5 Hz to 10 Hz. Improvement over the perfonnance of 
the embedded feedforward controller (figs. 35-41) is due to the incorporation of a feedforward 
path (the last two terms in eq. 15) and the ability of the explicit feedforward component (the last 
tenn in eq. 15) to track changes in the gust frequency. The closed-loop spectra (figs. 60-62) are 
reasonably flat, indicating that the controller is near optimal. The required control angles are 
quite large compared to that required for quelling the responses for constant-frequency gust 
excitations. As discussed earlier in the case of the embedded feedforward controller, this may be 
due to the system identification being compromised by the magnitude of the disturbance signal 
relative to the noise input on the control channels during identification. Again, it is suspected 
that if the controller were fully adaptive (a c , [f, S‘\ and / in equation 15 all adapt at each time 
step) performance would be much better and require much smaller control angles. 

Summary of Observations : Taken as a whole, the simulated closed-loop behavior of the HiLDA 
wing model during constant- and swept-frequency vertical gust excitation indicates that GPC- 
based controllers are able to significantly reduce wing bending and torsion moments and to do so 
rapidly after closing the loop. Tables 1 and 2 summarize the overall factors by which the SBI2, 
SBY, and STI2 wing root moments were reduced for each of the three gusts for the cases of 
embedded feedforward and explicit feedforward, respectively. For the case of embedded 
feedforward (table 1), the reductions range from a factor of four to as high as 667, depending on 
whether the excitation is from a constant- or swept-frequency gust and on the direction of the 
response. For the case of explicit feedforward (table 2), the reductions range from a factor of 
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four to as high as 1000, again depending on whether the excitation is from a constant- or swept- 
frequency gust and on the direction of the response. The control angles required to effect these 
reductions ranged from one to three degrees for fixed- frequency gust excitations and from 10 to 
20 degrees for swept-frequency gusts. Viewed as a whole, these results are decidedly positive in 
regard to establishing the credibility of GPC as a serious candidate approach for active GLA. 
Additional wind-tunnel testing of the HiLDA wing model is planned to provide additional 
validation of the method. 


Concluding Remarks 

The results of numerical simulations aimed at assessing and evaluating the efficacy of a 
predictive control method known as Generalized Predictive Control (GPC) for active gust load 
alleviation (GLA) using trailing- and leading-edge control surfaces have been summarized. GPC 
is a linear, time-invariant, multi-input/multi-output predictive control method that uses an 
AutoRegressive with Exogenous input (ARX) model to describe the input-output relationship of 
the system. The coefficient matrices of the ARX equation are determined using system 
identification techniques. The input-output equation is used to fonn a multi-step output 
prediction equation over a finite prediction horizon. The control to be imposed at the next time 
step is determined by minimizing the deviation of the predicted controlled plant outputs from the 
desired (or target) outputs, subject to a penalty on control effort. The method lends itself to an 
implementation that allows all computations to be done on-line and in real time. Following a 
qualitative discussion of the essential features of GPC, the equations underlying the method were 
presented and discussed, including system identification, calculation of control law matrices, and 
calculation of the commands applied to the control effectors. Both embedded and explicit 
feedforward paths for incorporation of disturbance effects were addressed. 

Representative results from two types of simulations were presented. The first used a 
three-degree-of-freedom mathematical model of a mass-spring-dashpot system subject to a 
variety of user-defined external disturbances. The second used open-loop data obtained in a 
wind-tunnel while subjecting a wing model to constant- frequency and swept-frequency 
sinusoidal vertical gust excitations. Closed-loop behavior of the model was simulated in post- 
test calculations. Results obtained from these simulations have been decidedly positive. In 
particular, the results of all the closed-loop simulations conducted for the wing model showed 
reductions in wing root moments ranging from factors of four to as much as 1000, depending on 
whether the excitation is from a constant- or swept-frequency gust and on the direction of the 
response. The control angles required to achieve these reductions ranged from one to three 
degrees for fixed-frequency gust excitations and from 10 to 20 degrees for swept-frequency 
gusts. The GPC-based controller was evaluated for both embedded feedforward and explicit 
feedforward modes of operation. The effectiveness of the explicit feedforward controller using 
an adaptive AutoRegressive (AR) model for predicting future disturbances was particularly 
notable, both in improving performance and in tracking the variable disturbance frequency. 

Based on the results of these studies, it appears that GPC is a strong candidate approach for 
active GLA and warrants further investigation in the subject area. Additional wind-tunnel testing 
of the HiLDA wing model is planned to provide additional validation of the method. 
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Appendix A 

Treatment of Measurable Coherent Disturbances 


If the disturbance is measurable and coherent, an AutoRegressive (AR) model of the 
disturbance giving the disturbance at a future time step k can be written in the form 

d(k) = rj x d(k-l)+rj 2 d(k-2)-\ \-rj n d(k-n d ) (Al) 

where the coefficient matrices //, are rj x rj, the disturbance vectors d are /v x /, r,i is the number 

of disturbances, and n ( i is the order of the AR model. The coefficient matrices //, in the model 

can be determined using established system identification techniques. The disturbance d(k+j ) at 
an arbitrary future time step k+j can be written as 

d(k+j) = rj l d(k+j-l)+rj 2 d(k+j-2) + --+i] nd d(k+j-n d ) (A2) 

Letting / range over the set of values 0 to h p - 1 one obtains the set of equations 


d(k) = rj x d(k-Y)+ij 2 d(k-2) J t vrj n d{k-n d ) 

d(k+l) = Jj l d(k)+Jj 2 d(k-l)d \-rj n d(k-n d + 1) 

d(k + 2 = 7) x d{k + 1)+ ?j 2 d(kd) -t f fj n d(k— n d +2) 

d(k + 3) — 7j x d{k +2 +T/ P d(k +1 ')-\-T/^d(k') ■+ f Tj n d{k — n d + 3) 


d(k + h p -1 ) = rj x d(k+h p -2)+ij 2 d(k + h p -3)d \-7] n ^d(k + h p —n d -1) 


The first equation in this set is the AR equation given in (Al) and succeeding equations are given 
using previously computed future disturbances (indicated in color) and measured past 
disturbances according to the pattern shown. The shifting pattern of the computations indicated 
in these equations is easily implemented in MATLAB. The results of these shifting calculations 
can be assembled into a matrix Df p relating future disturbances to past disturbances according to 
the equation 


{d f (k)} = [d fp ]{d p {k - p)} 

h r. x 1 h r , x pr , pr, x 1 

pa p a r a r d 

and the term D t df(k) in equation 10 can be replaced by D f Df p d p (k-p). 
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Batch Identification of Disturbance Models 


Comparisons of several imposed/measured disturbance time histories with those 
predicted using a non-adaptive AR model are summarized in figures A1-A5. The AR model was 
identified using a small block of “measured” disturbance data. These results are intended to 
provide an indication of the ability of AR models to adequately represent correlated disturbances 
for GPC-based active control systems. 

The imposed disturbance time history used for figures A1 and A2 is composed of four 
equal-length time segments, each containing 100 data points, for a total time period of 20 

seconds (4*100*0.05 = 20 sec, where At = 0.05 sec). Each segment has a different amplitude 

and frequency. The predictions in figure A1 were based on a model that was identified using 
only the first 50 data points (first 50*0.05 = 2.5 sec of data) while the predictions in figure A2 
were based on a model that used all 400 data points. A 6 th order (six-tenn) AR model {rid = 6) 
was assumed for both. These results suggest that a disturbance that has considerable variation in 
its amplitude and/or frequency can be represented rather well by a low-order model based on a 
single identification using a modest-length block of data if there is correlation in the disturbance. 
This appears to be so even if there are abrupt changes in the character of the disturbance (and 

hence in the coefficient matrices r/ ( in equation Al). 

Figures A3-A5 show comparisons of gust angles measured during the HiLDA test at a 
dynamic pressure of 51 lb/ft 2 with predictions obtained using a six-term AR model of the 
disturbance. The data were obtained at a sampling rate of 500 Hz but were decimated by a factor 

of ten to reduce the effective sampling rate to 50 Hz. Hence, the effective sampling time is At = 

0.02 sec. The decimated time history contained 1300 points. The first 500 data points (first 
500*0.02 = 10 sec) of this decimated time history were used to identify the coefficient matrices 

rii of the AR disturbance model. Figure A3 shows the measured versus predicted time histories 

for the case in which the gust frequency was held constant near the first vertical bending mode 
(about 2 Hz). Similar results for a constant-frequency excitation near the second vertical 
bending mode (about 8 Hz) are given in figure A4. Figure A5 is for the case in which the gust 
frequency was swept from 0 to 10 Hz in 30 seconds. In spite of the fact that the gusts are not 
pure sinusoids and are contaminated by higher hannonics and noise, there appears to be 
sufficient coherence in the signals to allow use of non-adaptive methods to identify disturbance 
models that produce reasonable results, particularly at the lower frequencies. 

Adaptive Identification of Disturbance Models 

Sometimes it is necessary to update the model of a system or a disturbance on-line and in 
real time while a system is in operation. Methods using models that are adjusted on-line using 
measured input/output data up to the current time are called adaptive. To be adaptive in real time 
requires that the measurements obtained at any time point can be processed during one sampling 
period. In such instances, recourse must be made to recursive identification techniques, which 
process data sequentially as they become available. A variety of recursive estimation methods 
are described in the literature (see, for example, refs. 11 and 21). Among these, the classical 
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recursive least squares (RLS) method is probably the most popular and simplest. For models 
whose assumed order is relatively small (such as the disturbance models of interest here), it is the 
method of choice. However, there are several fast versions of the method available that would 
be better suited to large models. 

Adaptive identification of a disturbance model would update the r ( i x r ( / coefficient 
matrices 77 , in equation A1 at each time step. This updated model would then be used at each 

time step to compute the vector df(k ) in equation A4 and then df(k) used in equations 10, 11, 12, 
14, and 15 as appropriate. 

Although the results of figures A1-A5 indicate that a non-adaptive identification scheme 
may be adequate for a wide range of disturbance time histories, there could be times when an 
adaptive scheme may prove advantageous or necessary. Figures A6-A9 show comparisons of 
several imposed and measured disturbance time histories with those predicted using an AR 
model of the disturbance that is updated at each time step using the classical recursive least 
squares technique. Figures A 6 and A7 show results using 2-term and 6 -term AR models, 
respectively. The imposed disturbance time history used is the same one that was employed in 
figures A1 and A2. Increasing the order of the AR model gives a predicted disturbance that is in 
better agreement with the imposed disturbance, except when transitioning from one type of 
disturbance to another where the differences increase away from transition points. 

Figures A 8 and A9 show results for the case of random disturbances using 6 -term and 
16-term AR models, respectively. No agreement is evident (or should be expected) because 
there is no correlation in random signals. 
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Figure Al.- Imposed versus predicted disturbance time histories using a 6-tenn AR model 
identified using the first 50 data points (first 2.5 sec of data). 
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Figure A2.- Imposed versus predicted disturbance time histories using a 6-tenn AR model 

identified using all 400 data points. 
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Figure A3.- Measured versus predicted gust time histories using a 6-term AR model identified 
using the first 500 data points (out of 1300) from HiLDA wing test. Airstream was excited 

at 2 Hz. 
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Figure A4.- Measured versus predicted gust time histories using a 6-term AR model identified 
using the first 500 data points (out of 1300) from HiLDA wing test. Airstream was excited 


at 8 Hz. 
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Figure A5.- Measured versus predicted gust time histories using a 6-tenn AR model obtained 
using the first 500 data points (out of 1300) from HiLDA wing test. Airstream was excited by 

sweeping gust frequency from 0.5 to 10 Hz. 
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Figure A6.- Imposed versus predicted disturbance time histories using a 2-tenn AR model 

updated at each time step using RLS. 
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Figure A7.- Imposed versus predicted disturbance time histories using a 6-tenn AR model 

updated at each time step using RLS. 
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Figure A8.- Imposed versus predicted disturbance time histories using a 6-tenn AR model 

updated at each time step using RLS. 
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Time, sec 

Figure A9.- Imposed versus predicted disturbance time histories using a 16-term AR model 

updated at each time step using RLS. 
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Appendix B 


State-Space Equations for 3-DOF Mass-Spring-Dashpot System 


The second-order matrix equations of motion for the 3-degree-of-freedom mass-spring- 
dashpot system shown in figure 2 are given by 


m l 0 0 [ Jc, I c, + c 2 -c. 

0 m 2 0 u, + —c 2 c 2 + 

0 0 m 3 x 3 J 0 -c 

ud { 
ud 2 
ltd, 

UC\ 

uc 2 
uc 3 

where the external forces acting on the system are divided into disturbance forces and control 
forces u C j for convenience in controls simulation work. The matrix \Bi\ is an input influence 
matrix indicating the locations and types of force inputs. Equation B1 can be written compactly 
as 

[M]{x] + [C]{i} + [£]{*} = {/} = [ S 2 ]{» } (B2) 

The corresponding first-order equations describing the continuous-time state-space model 
for the system are given by the state equation 

{VH-ilW + [£]{«} (B3) 

and the output equation 

W = [Cj{X} + [Z) c ]{n} (B4) 


"1 0 0 1 0 0 " 
= 010010 
0 0 1 0 0 1 



where the state matrix [A c ], state matrix [Zf], and state vector {X} are given by 




[o] [/] 

-[at'x] -[at'c] ’ 


[«.]= 


[0] 

\m-%y 
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The outputs (measured physical quantities) are related to the state variables and the inputs by the 
output equation which can be written as 

M = [C„ C„]{^)} + [Cj{x} (B6) 

where Q, C v , and C a are 3><3 output influence matrices for displacement, velocity, and 
acceleration, respectively, and have the form 



These matrices describe the relationship between the vectors {x},{x},{x} and the measurement 

vector {y}. If the measurement unit is a physical unit (e.g., displacement, velocity, or 
acceleration) the corresponding tenn is unity; otherwise there is a conversion factor from 
physical units to the measurement units. 

Solving equation B2 for {x} and substituting into equation B6 gives 

bHC„ Cj{«} + [Cj[-[M-.C] W -[M-,] { q + [M-., ; ]W] (B8a) 

or 

{y} = [C d -C„M- l K C,,-C a M- | c]||*|| + [c,,M- | i5 ! ]{ u } (B8b) 

or 

{y}=[C ll -C ll M-'K C,-C a M- | c]{X} + [c„M* 1 S 2 ]{»} (B8c) 

which can be written compactly as 

M=[cJM + [£JM (B9) 

where \C c ] and [ D ] are the output and direct transmission matrices and are given by 
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(BIO) 


[C c ] = [C d -C a M-'K C V -C„M-'C 




The direct transmission matrix [/) ] is zero if accelerations are not included in the output 

measurements {y } (in which case [C a ] = [0] ), or if accelerations are measured but there are no 
external forces (disturbance and control) acting at those locations (in which case some of the 
diagonal elements of I u d and I uc will be zero). 

Special cases of the general equations given above are easily obtained. For example, the 
number, type, and combination of input and output quantities can be varied by deleting 
appropriate rows of the output equation and/or zeroing appropriate diagonal elements in I u d, I uc , 
Cdi f .. and C a . 
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Table 1.- Reduction Factors for HiLDA Wing Model Moment Responses 

Using Embedded Feedforward 


Gust Frequency 

SBI2 

SBY 

STI2 

2 Hz 

667 

10 

34 

8 Hz 

75 

33 

30 

Swept from 0.5 - 10 Hz 

31 

7 

4 


Table 2.- Reduction Factors for HiLDA Wing Model Moment Responses 

Using Explicit Feedforward 


Gust Frequency 

SBI2 

SBY 

STI2 

2 Hz 

1000 

30 

170 

8 Hz 

500 

100 

60 

Swept from 0.5 - 10 Hz 

625 

40 

4 
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Figure 1.- Block diagram of GPC identification and control procedure. 
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Figure 2.- Three-degree-of-freedom mass-spring-dashpot system. 
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Figure 3.- Steady-state sinusoidal disturbances imposed on 3-DOF math model 

at masses 1 and 2 {udi and udi). 



(a) Piecewise varying frequency and amplitude 



(b) Continuously varying frequency and amplitude 
Figure 4 - Non-steady-state disturbances imposed on 3-DOF math model at mass 2 (uch). 
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Figure 5.- Feedback only for 3 DOF system with steady-state disturbances (Control: uc 3 ; 

Output: x l ; Disturbance: udj and udj). 
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Figure 6.- Feedback with embedded feedforward for 3 DOF system with steady-state 
disturbances (Control: wc Output: x l ; Disturbance: udi and udi). 
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Figure 7.- Feedback with explicit feedforward of both steady-state disturbances for 3 DOF 
system (Control: ucf, Output: x l ; Disturbance: udi and udj). 
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Figure 8.- Feedback with explicit feedforward of both steady-state disturbances for 3 DOF 
system (Control: uci and uc f. Output: x l ; Disturbance: udi and udi). 
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Figure 9.- Feedback with explicit feedforward of both steady-state disturbances for 3 DOF 
system (Control: uci and UC 3 ; Output: x l and x 2 ; Disturbance: udi and udi) 
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Figure 10.- Feedback with embedded feedforward of disturbance with piecewise varying 
amplitude and frequency (Control: wc?; Output: x l ; Disturbance: udi). 
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Figure 11.- Feedback with explicit feedforward of disturbance with piecewise varying amplitude 
and frequency (Control: ucf, Output: x l ; Disturbance: udi). 
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Figure 12.- Feedback with embedded feedforward of disturbance with piecewise varying 
amplitude and frequency (Control: uci and UC 3 ; Output: x l ; Disturbance: udi). 
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Figure 13.- Feedback with explicit feedforward of disturbance with piecewise varying amplitude 
and frequency (Control: uci and UC 3 ; Output: x l ; Disturbance: udi). 
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Figure 14.- Feedback with explicit feedforward of disturbance with piecewise varying amplitude 
and frequency (Control: uci, UC 2 , and ucj; Output: x,; Disturbance: udi). 
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Figure 15.- Feedback with embedded feedforward of disturbance with continuously varying 
frequency and amplitude (Control: wc?; Output: .y; Disturbance: udi). 
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Figure 16.- Feedback with explicit feedforward of disturbance with continuously varying 
frequency and amplitude acting at mass 2 (Control: UC 3 ; Output: x l ; Disturbance ad 2 ). 


53 



Time, sec 



Figure 17.- Feedback with explicit feedforward of disturbance with continuously varying 
frequency and amplitude (Control: uci and uci. Output: x l ; Disturbance: udi). 
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Figure 18.- HiLDA wing model in Langley Transonic Dynamics Tunnel. View looking 

upstream. 


55 



Figure 19.- Biplane vanes of TDT airstream oscillator system. View looking downstream toward 

model in test section. 
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Figure 20.- Schematic of HiLDA model wing planform showing arrangement of control surfaces 
and locations of strain-gage response sensors used for simulations. 
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Figure 21.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 
gust excitation at 2 Hz. Embedded feedforward. 
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Figure 22.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) during 

gust excitation at 2 Hz. Embedded feedforward. 
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Figure 23.- Time histories of open- and closed-loop torsional moment (STI2) during 
gust excitation at 2 Hz. Embedded feedforward. 
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Figure 24.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 
during gust excitation at 2 Hz. Embedded feedforward. 
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Figure 25.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 
(SBY) during gust excitation at 2 Hz. Embedded feedforward. 
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Figure 26.- Open- and closed-loop power spectral densities for torsional moment (STI2) 
during gust excitation at 2 Hz. Embedded feedforward. 
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Figure 27.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust excitation 

at 2 Hz. Embedded feedforward. 
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Figure 28.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 
gust excitation at 8 Hz. Embedded feedforward. 
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Figure 29.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) 
during gust excitation at 8 Hz. Embedded feedforward. 
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Figure 30.- Time histories of open- and closed-loop torsional moment (STI2) during 
gust excitation at 8 Hz. Embedded feedforward. 
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Figure 3 1 Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 
during gust excitation at 8 Hz. Embedded feedforward. 
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Figure 32.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 
(SBY) during gust excitation at 8 Hz. Embedded feedforward. 
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Figure 33.- Open- and closed-loop power spectral densities for torsional moment (STI2) 
during gust excitation at 8 Hz. Embedded feedforward. 
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Figure 34.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust 
excitation at 8 Hz. Embedded feedforward. 
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Figure 35.- Time histories of open- and closed-loop vertical bending moment (SBI2) while 
sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 


72 



Figure 36.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) while 
sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
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Figure 37.- Time histories of open- and closed-loop torsional moment (STI2) while sweeping 
gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
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Figure 38.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 
while sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
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Figure 39.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 
(SBY) while sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
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Figure 40.- 


Open- and closed-loop power spectral densities for torsional moment (STI2) while 
sweeping gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
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Figure 41.- Time histories of closed-loop control inputs (CTL1-CTL5) while sweeping 
gust frequency from 0.5 to 10 Hz. Embedded feedforward. 
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Figure 42.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 

gust excitation at 2 Hz. Explicit feedforward. 


79 



Figure 43.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) during 

gust excitation at 2 Hz. Explicit feedforward. 
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Figure 44.- Time histories of open- and closed-loop torsional moment (STI2) during 
gust excitation at 2 Hz. Explicit feedforward. 
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Figure 45.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 
during gust excitation at 2 Hz. Explicit feedforward. 


82 


SBY, (in-lb) /Hz 



Figure 46.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 
(SBY) during gust excitation at 2 Hz. Explicit feedforward. 
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Figure 47.- Open- and closed-loop power spectral densities for torsional moment (STI2) 
during gust excitation at 2 Hz. Explicit feedforward. 
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Figure 48.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust 
excitation at 2 Hz. Explicit feedforward. 
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Figure 49.- Time histories of open- and closed-loop vertical bending moment (SBI2) during 

gust excitation at 8 Hz. Explicit feedforward. 
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Figure 50.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) during 

gust excitation at 8 Hz. Explicit feedforward. 
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Figure 51.- Time histories of open- and closed-loop torsional moment (STI2) during 
gust excitation at 8 Hz. Explicit feedforward. 
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Figure 52.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 
during gust excitation at 8 Hz. Explicit feedforward. 
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Figure 53.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 
(SBY) during gust excitation at 8 Hz. Explicit feedforward. 
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Figure 54.- Open- and closed-loop power spectral densities for torsional moment (STI2) 
during gust excitation at 8 Hz. Explicit feedforward. 
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Figure 55.- Time histories of closed-loop control inputs (CTL1-CTL5) during gust 
excitation at 8 Hz. Explicit feedforward. 
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Figure 56.- Time histories of open- and closed-loop vertical bending moment (SBI2) while 
sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
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Figure 57.- Time histories of open- and closed-loop fore-and-aft bending moment (SBY) 
while sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
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Figure 58.- Time histories of open- and closed-loop torsional moment (STI2) while sweeping 
gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
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Figure 59.- Open- and closed-loop power spectral densities for vertical bending moment (SBI2) 
while sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
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Figure 60.- Open- and closed-loop power spectral densities for fore-and-aft bending moment 
(SBY) while sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
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Figure 61.- Open- and closed-loop power spectral densities for torsional moment (STI2) while 
sweeping gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
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Figure 62.- Time histories of closed-loop control inputs (CTL1-CTL5) while sweeping 
gust frequency from 0.5 to 10 Hz. Explicit feedforward. 
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